
#define MIN(A,B) ( (A) < (B)? (A) : (B)  )
#define MAX(A,B) ( (A) > (B)? (A) : (B)  )
#define M2V(i, j, LDA) ((j) * (LDA) + (i))

#define EPSILON 1.0e-16

/* Maximum number of rows or cols allowed in print_matrix */
#define PRINT_LIMIT 0


/**** remove me ****/
typedef enum{
    CUBLAS_STATUS_SUCCESS         =0,
    CUBLAS_STATUS_NOT_INITIALIZED =1,
    CUBLAS_STATUS_ALLOC_FAILED    =3,
    CUBLAS_STATUS_INVALID_VALUE   =7,
    CUBLAS_STATUS_ARCH_MISMATCH   =8,
    CUBLAS_STATUS_MAPPING_ERROR   =11,
    CUBLAS_STATUS_EXECUTION_FAILED=13,
    CUBLAS_STATUS_INTERNAL_ERROR  =14
} cublasStatus_t;

struct cublasContext;
typedef struct cublasContext *cublasHandle_t;

#define CUBLASAPI 

typedef enum {
    CUBLAS_FILL_MODE_LOWER=0, 
    CUBLAS_FILL_MODE_UPPER=1
} cublasFillMode_t;

typedef enum {
    CUBLAS_DIAG_NON_UNIT=0, 
    CUBLAS_DIAG_UNIT=1
} cublasDiagType_t; 

typedef enum {
    CUBLAS_SIDE_LEFT =0, 
    CUBLAS_SIDE_RIGHT=1
} cublasSideMode_t; 


typedef enum {
    CUBLAS_OP_N=0,  
    CUBLAS_OP_T=1,  
    CUBLAS_OP_C=2  
} cublasOperation_t;


typedef enum { 
    CUBLAS_POINTER_MODE_HOST   = 0,  
    CUBLAS_POINTER_MODE_DEVICE = 1        
} cublasPointerMode_t;



   /* acc_init( devicetype )
   Initialized the runtime system and sets the accelerator device
   type to use for this host thread.
   Usage example 
   acc_init(acc_device_type);
   */
 /*  typedef  enum  {
      acc_device_none,
      acc_device_default,
      acc_device_host,
      acc_device_not_host,
      acc_device_nvidia,
      acc_device_generic_opencl
   } acc_device_t;*/



 
/* Cholesky factorization (outer product version)
   ===============================================

     Input: R  , matrix to be factorized
            N  , order of the matrix
            LDA

     Output: R ,  Upper matrix of the factorization
*/
int acc_cholesky_outer(double * R, int n, int lda) {
	int i;
	int correct = 1;
    // This is not the optimal OpenACC implementation. Probably it would be better to
    //   split this into smaller kernels and avoid using for inside kernel. 
   // **** Discuss with ivan #pragma acc kernels loop pcopyout(R[n*n]) private(i) firstprivate(lda)
   // #pragma acc kernels loop private(i) firstprivate(lda)
	for (i = 0; i < n; i++) {
		{
            int k;
			R[M2V(i, i, lda)] = sqrt(R[M2V(i, i, lda)]);
            
			for (k = (i + 1); k < n; k++) {
				R[M2V(i, k, lda)] = R[M2V(i, k, lda)] / R[M2V(i, i, lda)];
			}
			if (R[M2V(i, i, lda)] <= 0.0f) {
				correct = 0;
			}
			if (correct) {
				for (k = i + 1; k < n ; k++) {
					for (int j = i + 1; j < n; j++) {
						R[M2V(k, j, lda)] = R[M2V(k, j, lda)] - R[M2V(i, j, lda)] * R[M2V(i, k, lda)];
					}
				}
			}
		}
	}

	return correct;
}

/* Cholesky factorization by blocks


     Input: R  , matrix to be factorized
            N  , order of the matrix
            block_size 

     Output: R ,  Upper matrix of the factorization

*/
int acc_cholesky_blocks(double * R, int n, int block_size) {
	int i, j, k;
	int lda = n;   // Just for now
	int n_block = (n / block_size) + (n % block_size?1:0);	
//    acc_init(acc_device_nvidia);
    /* Initialize runtime */
    // This should create an internal scope
    // Cublas initialization
    cublasStatus_t status;
    cublasHandle_t handle;
    status = cublasCreate(&handle);
    // Frangollo gives access to its internal cuda stream
    //  This functionality is available also for OpenCL, 
    //   but there is not an standard call for this
    // This is not supported by any other compiler, again AFAIK
//    cudaStream_t = FRG__getCmdQueue(); 
    // AFAIK, If using an stream, cublas calls are asyncronous 
//    status = cublasSetStream(handle, stream);
     
	{
	for (k = 0; k < n_block; k++) {
        int ii, jj;
		int bs_k = MIN(block_size, n - k * block_size);
        double * cholBlockDiscont = R + M2V(k * block_size,k * block_size,lda);
        double * cholBlock = (double *) malloc(block_size * block_size * sizeof(double));
        for (ii = 0; ii < bs_k; ii++)
            for (jj = 0; jj < bs_k; jj++) {
                 cholBlock[ii * bs_k + jj] = cholBlockDiscont[ii * lda + jj];
             }
                   /* esta rutina tiene un kernels, miralo */
                if (!acc_cholesky_outer(cholBlock, bs_k, bs_k)) {
                     printf("*** Cholesky for submatrix [%d][%d] failed \n", k, k);
                }
     #pragma acc data copy(cholBlock[block_size * block_size])
        {
                double * B = (double *) malloc(block_size * block_size * n_block * sizeof(double));
                double * dtrsm_init = R + M2V(k * block_size,(k+1) * block_size,lda);
                int ii, jj;
                for (ii = 0; ii < (n - ((k+1) * block_size)); ii++) 
                    for (jj = 0; jj < bs_k; jj++)
                        B[ii * block_size + jj] = *(dtrsm_init + M2V(ii,jj,lda));

            #pragma acc data copy(B[bs_k*n_block*block_size]) 
                   {
                /** Copy the entire n_block*bs **/
                for (j = k + 1; j < n_block; j++) {
                    double ALPHA = 1.0;
                    int bs_j = MIN(block_size, n - j * block_size);
                  #pragma acc host_data use_device(cholBlock, B)  
                        cublasDtrsm(handle, CUBLAS_FILL_MODE_LOWER, 
                                    CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_T, 
                                    CUBLAS_OP_N, bs_k, bs_j, &ALPHA, 
                                    (double *) cholBlock, bs_k, 
                                    (double *) (B + (j * block_size)), block_size ) ;
                        
                        }
                
                   } /* copy B s */
               for (ii = 0; ii < (n - ((k+1) * block_size)); ii++) 
                            for (jj = 0; jj < bs_k; jj++)
                                (*(dtrsm_init + M2V(ii,jj,lda))) = B[ii * block_size + jj];

                    free(B);


    /**** Prepare dgemm **/
                for (i = k + 1; i < n_block ; i++) {
                                double ALPHA = -1.0;
                                double BETA = 1.0;

                      int bs_i = MIN(block_size, n - i * block_size);
                      double * fA = R + M2V(k * block_size, i * block_size, lda);
                      double * A = (double * ) malloc(block_size * block_size * sizeof(double));
                        for (ii = 0; ii < block_size; ii++)
                                for (jj = 0; jj < block_size; jj++) {
                                        A[ii * block_size + jj] = fA[ii * lda + jj];
                                }
                    #pragma acc data copyin(A[block_size*block_size])
                    {
                            for (j = i; j < n_block; j++) {
                                // ----- Scalar R[M2V(i, j, n)] = R[M2V(i, j, n)] - R[M2V(k, j, n)] * R[M2V(k, i, n)];
                                // ----- Block Rij = Rij - Rki^T Rkj
                                // Blas:
                                //  1. Rij = Rij - Rki^T * Rkj
                                int bs_j = MIN(block_size, n - j * block_size);
                    //			dgemm_("T", "N", &bs_i, &bs_j, &bs_k, &ALPHA, R + M2V(k * block_size, i * block_size, lda), &lda, R + M2V(k * block_size, j * block_size, lda), &lda, &BETA, R + M2V(i * block_size, j * block_size, lda),  &lda);
                                double * fB = R + M2V(k * block_size, j * block_size, lda);
                                double * fC = R + M2V(i * block_size, j * block_size, lda);
                                double * B = (double * ) malloc(block_size * block_size * sizeof(double));
                                double * C = (double * ) malloc(block_size * block_size * sizeof(double));
                            for (ii = 0; ii < block_size; ii++)
                                for (jj = 0; jj < block_size; jj++) {
                                     B[ii * block_size + jj] = fB[ii * lda + jj];
                                     C[ii * block_size + jj] = fC[ii * lda + jj];
                                     }

                       #pragma acc data copyin(B[block_size*block_size]) copy(C[block_size*block_size]) 
                           {
                             
                               #pragma acc host_data use_device(A, B, C)  
                                status = cublasDgemm(handle, CUBLAS_OP_T, CUBLAS_OP_N, 
                                                    bs_i, bs_j, bs_k, &ALPHA, 
                                                    (double *)A, block_size, (double *)B, 
                                                    block_size, &BETA, (double *)C ,  block_size);
                            }
                              for (ii = 0; ii < block_size; ii++)
                                for (jj = 0; jj < block_size; jj++) {
                                     fC[ii * lda + jj] = C[ii * block_size + jj];
                                     }

                            free(B);free(C);
                            } /* block mxm */ 
#if 0
                                /* cublasDsyrk_("U", "T", &bs_i, &bs_k, &ALPHA, 
                                    R + M2V(k * block_size, i * block_size, lda), &lda, &BETA, 
				                	R + M2V(i * block_size, i * block_size, lda), &lda); */
                      double * fii = R + M2V(i * block_size, i * block_size, lda);
                      double * Rii = (double * ) malloc(block_size * block_size * sizeof(double));
                        for (ii = 0; ii < block_size; ii++)
                                for (jj = 0; jj < block_size; jj++) {
                                        Rii[ii * block_size + jj] = fii[ii * lda + jj];
                                }
                        #pragma acc data copy(Rii[block_size*block_size])
                        {
                               #pragma acc host_data use_device(A, Rii)  
                               cublasDsyrk(handle, CUBLAS_FILL_MODE_UPPER, CUBLAS_OP_T, 
                                           bs_i, bs_k, &ALPHA, 
                                        (double *) A, block_size, &BETA, 
	    			                	(double *) Rii, block_size); 
                        } /* data Rii */
                        for (ii = 0; ii < block_size; ii++)
                                for (jj = 0; jj < block_size; jj++) {
                                        fii[ii * lda + jj] = Rii[ii * block_size + jj];
                                }
                        free(Rii);
#endif
                        } /* data A */
                        free(A);
              
                 }  /* tailing submatrix loop end */
                 
		            } /* data */
        for (ii = 0; ii < bs_k; ii++)
            for (jj = 0; jj < bs_k; jj++) {
                 cholBlockDiscont[ii * lda + jj] = cholBlock[ii * bs_k + jj];
             }
		}
	}
    FRG__printStats();
    // Destroy here all cublas stuff
    // This should destroy runtime context and clear scope
    //    cublasShutdown();
    cublasDestroy(handle);
	return 1;
}

